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ABSTRACT 

We investigate orbital resonances expected to arise when a system of two planets, 
with masses in the range 1-4 M§ , undergoes convergent migration while embedded in 
a section of gaseous disc where the flow is laminar. We consider surface densities cor- 
responding to 0.5 — 4 times that expected for a minimum mass solar nebula at 5.2 AU. 
For the above mass range the planets undergo type I migration. Using hydrodynamic 
simulations we find that when the configuration is such that convergent migration 
occurs the planets can become locked in a first order commensurability for which the 
period ratio is (p + l)/p with p being an integer and migrate together maintaining it 

• for many orbits. Slow convergent migration results in commensurabilities with small 
C""'- p such as 1 or 2. Instead, when the convergent migration is relatively rapid as tends 

to occur for disparate masses, higher p commensurabilities are realized such as 4:3, 
5:4, 7:6 and 8:7. However, in these cases the dynamics is found to have a stochastic 
character with some commensurabilities showing long term instability with the con- 

• sequence that several can be visited during the course of a simulation. Furthermore 
Q^' the successful attainment of commensurabilities is also a sensitive function of initial 

conditions. When the convergent migration is slower, such as occurs in the equal mass 
£^ ■ case, lower p commensurabilities such as 3:2 are obtained which show much greater 

^ | stability. 

Resonant capture leads to a rise in eccentricities that can be predicted using a simple 
analytic model, that assumes the resonance is isolated, constructed in this paper. We 
find that, once the commensurability has been established, the system with an 8:7 
. commensurability is fully consistent with this prediction. 

' We find that very similar behaviour is found when the systems are modeled using an N 

body code with simple prescriptions for the disc planet interaction. Comparisons with 
the hydrodynamic simulations indicate reasonably good agreement with predictions 
for these prescriptions obtained using the existing semi-analytic theories of type I 
migration. 

We have run our hydrodynamic simulations for up to 10 — 10 5 orbits of the inner 
planet. Longer times could only be followed in the simpler N-body approach. Using 
that, we found that on the one hand an 8:7 resonance established in a hydrodynamic 
simulation could be maintained for more than 10 6 orbits. On the other hand other 
similar cases show instability leading to another resonance and ultimately a close 
scattering. 

There is already one known example of a system with nearly equal masses in the 
several Earth mass range, namely the two pulsar planets in PSR B1257+12 which are 
intriguingly, in view of the results obtained here, close to a 3:2 commensurability. This 
will be considered in a future publication. 

Future detection of other systems with masses in the Earth mass range that display 
orbital commensurabilities will give useful information on the role and nature of orbital 
migration in planet formation. 
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1 INTRODUCTION 

The increasing number of extrasolar multi-planet systems, 
their diversity, and dynamical complexities provide a strong 
motivation to study the evolution and stability of such sys- 
tems. One of the important features connected with plane- 
tary system evolution is the occurrence of mean motion reso- 
nances, which may relate to conditions at the time of or just 
after the process of formation. There are some well known 
examples of systems of celestial bodies exhibiting mean mo- 
tion resonances both within our Solar System (Neptune and 
Pluto; Io, Ganymede and Europa (see eg. Goldreich, 1965) ) 
and outside such as Gliese 876 (Marcy et al., 2001), HD 
82943 (Mayor et al, 2004) and 55 Cancri (Mc Arthur et al, 
2004). As the latter examples involve planets with masses 
in the Jovian mass range, most of the investigations appro- 
priate to extrasolar planets have focused on giant planets. 
However, it is likely that planetary systems around other 
stars may harbor planets with masses in the Earth mass 
range as well. These should be revealed by future space- 
based missions, such as Darwin, COROT, Kepler, SIM and 
TPF. 

Meanwhile, it is important to establish the main 
features of the evolution of low mass planets embed- 
ded in a gaseous disc and in particular to determine 
the types of resonant configurations that might arise 
when a pair of such planets evolves together. The disc 
planet interaction naturally prod uces orbital migration 
through the action of tidal tor ques dGoldreich fe Tremaind 
Il980l : iLin fe Papaloizoul Il986l) which in turn may lead 
to an orbital resonance in a many planet system eg . 
dSnellgrove. Papaloizou fc Nelson! l200ll : iLee fc Pealel 120021) . 
For low mass planets the disc undergoes small linear pertur- 
bations that induce density waves that propagate away from 
the planet. The angular momentum these waves transport 
away results in r apid orbital migration called type I migra- 
tion dWardll997l) . In this type of migration, when the disc is 
laminar and inviscid, the planet is embedded and the surface 
density profile of the disc remains approximately unchanged. 
The rate of migration is proportional to mass of the planet 
and the timescale of inward migration on a circular orbit 
can be estimate d for a disc with constant surface density to 
be given by (see lTanaka. Takeuchi fc WardI (120021) ) 



W„ 



(i) 



Here M* is the mass of the central star, m p i anet is the mass 
of the planet orbiting at distance r — r p , E p is the disc 
surface density at r = r p , c and Q p are the the local sound 
speed and angular velocity, Q, at r = r p respectively. The 
numerical coefficient W m = 0.3704. 

It is important to note that type I migration appropriate to 
a laminar disc may lead to short migration times in stan- 
dard model discs, that may th reaten the survival of pro- 
toplanetary cores dWardl dl997D ). However, in a disc with 
turbulence driven by the magnetorotational instability, the 
migration may be sto c hastic and accordingly less effective 
dNelson fc Papaloizou! J2004)). Nonetheless there is consid- 
erable uncertainty as to the extent of turbulent regions in 
the disc resulting from uncertainties in the degree of ioniza- 
tion (eg. Fromang, Terquem & Balbus, 2002) so that type 
I migration appropriate to a laminar disc may operate in 



some regions. We also note that because it is inversely pro- 
portional to the disc surface density, the migration time be- 
comes long in low surface density regions. Accordingly, for 
this first study of resonant interactions of planets in the 
mass range 1 — 3OM0 embedded in a gaseous disc, we shall 
consider only migration induced by a laminar disc. 
The density waves excited by a low mass planet with 
small eccentricity also lead to orbital circulariz ation (eg. 
lArtvmowiczl dl993l) : lPapaloizou fc L arwoodJ (|2000j)) at a rate 
that can be estimated to be given bv lTanaka fc WardI tQQ4\ 



tc = 



W c \ r p Q, p 



(2) 



Here, the numerical coefficient W c = 0.289. 

It is expected from equation Q that two planets with 
different masses will migrate at different rates. This has 
the consequence that their period ratio will evolve with 
time and may accordingly attain and, in the situation 
where the migration is such that the orbits converge, sub- 
sequently become locked in a mean motion resonance (e g. 
iNelson fc Panaloizoul (|2002l) : lKlev. Peitz fc Brvdenl l|2004|) ). 
In the simplest case of nearly circular and coplanar orbits 
the strongest resonances are the first-order resonances which 
occur at locations where the ratio of the two orbital periods 
can be expressed as the ratio of two consecutive integers, 
(p + l)/p, with p being an integer. As p increases, the two 
orbits approach each other and the strength of the reso- 
nance increases. In addition, the distance between succes- 
sive resonances decreases as p increases. The combination 
of these effects ultimately causes successive resonances to 
overlap and so, in the absence of gas, leads to the onset of 
chaotic motion. 

Resonance overlap occurs when the difference of the semi- 
major axes of the two planets is below a limit with half-width 
given, in the case of two equal mass planets, by iGladmanl 



Aa ^ 2_ ^ 2 ( m P lane 
a ~ 3p ~ V M* 



2/7 



(3) 



with a and m p ; anet being the mass and semi-major axis of 
either one of them respectively. Thus for a system consisting 
a two equal four Earth mass planets orbiting a central solar 
mass we expect resonance overlap for p > 8. Conversely we 
might expect isolated resonances in which systems of planets 
can be locked and migrate together if p < 8. 
However, note that the above discussion does not incor- 
porate the torques producing convergent migration or ec- 
centricity damping and thus may not give a complete 
account of the forms of chaos th at might be expected. 
iKarv. Lissauer fc Greenzweid dl993l) have discussed the case 
of small particles migrating towards a much more mas- 
sive planet and indeed conclude that chaotic behaviour is 
more extensive in the non conservative case. This occurs 
because the higher p resonances can be unstable thus pre- 
venting long term trapping. The instability arises because, 
if locked in resonance, the orbit of a planet has some ec- 
centricity. Making a slight perturbation to the orientation 
of the orbit produces an impulsive change to the semi- 
major axis at the next conjunction. This in turn affects 
the phase of the next conjunction and the next impul- 
sive change to the semi-major axis. For resonances of high 
enough p this sequence results in instability and chaotic be- 
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haviour. iKarv. Lissauer fc Greenzweid il993T> give an esti- 
mate of when this form of instability occurs as 



p > 0.667 



{ 8nem p i anet \ 1/5 
V 3M» J 



(4) 



Here m p i a net/M, refers to the mass ratio for the most mas- 
sive planet, and e is the orbital eccentricity for the planet 
with the smallest mass. Taking values appropriate to our 
simulations of e ~ 0.04, and m p i anet /M t ~ 10 -5 , gives 
p > 8. This is similar to the non dissipative estimate. 
When this form of instability operates, long term stable 
trapping in commensurabilities is not possible. However, a 
system can remain in one for a long time before moving into 
another higher p commensurability. Furthermore detailed 
outcomes are very sensitive to input parameters, which is 
characteristic of chaotic motion. Slight changes can alter the 
sequence of commensurabilities a system resides in making 
the issue of their attainment acquire a probabilistic char- 
acter. Both our hydrodynamic simulations and N body cal- 
culations show evidence of this behaviour. Although there 
may be islands of apparent stability, a system in this regime 
may ultimately undergo a scattering and exchange of the 
orbits of the two planets so that we should focus on sta- 
ble lower p commensurabilities as being physically possible 
planetary configurations. The arguments given above sug- 
gest that these must have p < 8. Our calculations indicate 
that the limit is even smaller. 

Hydrodynamic simulations of disc planet interactions, in 
which the discs are modelled as flat two dimensional objects 
with laminar flow governed by the Navier Stokes equations 
and which incorporate a migrating two giant planet sys- 
tem which e volves into a 2:1 commensurability have been 
perfor med bv lKlevI ll2000l ). ISnellgrove. Papaloizou fc Nelson! 
(2001) and successfully applied to the GJ876 system. In 
additi on to performing additional simulations, IPapaloizoul 
( 2003) has developed an analytic model describing two plan- 
ets migrating in resonance with arbitrary eccentricity. In this 
model the eccentricities are determined as a result of the 
balance between migration and orbital circularization. How- 
ever, disc tides were considered to act on the outer planet 
only. Capture of giant planets into re sonance has also been 
recent ly investigated numerically by iKlev. Peitz fc Brvdenl 
J2004D . 

In this paper we extend the types of two planet systems 
considered to include those with orbital resonances formed 
when the masses of the planets are 1 — 30M^ . Then for discs 
with aspect ratio H/r — c/(rQ) ~ 0.05 as in standard pro- 
tostellar disc models, the orbital migration will be induced 
by type I migration. 

In Section 2 we describe initial set up for our calculations 
giving the properties of the system in which two planets 
interact with each other and a gaseous disc. 
In Section 3 we present the results of two dimensional nu- 
merical calculations showing how convergent differential mi- 
gration of two planets can lead to resonance trapping fol- 
lowed by migration in which the resonance is maintained 
for the duration of the simulation which is characteristically 
around 2 x 10 4 orbits. By considering simulations with a 
range of planet masses and initial conditions, we are able to 
explore a variety of commensurabilities and show the sensi- 
tivity of the attainment of higher p commensurabilities to 



the input parameters, this being indicative of stochastic be- 
haviour. 

In Section 4 we compare our results with simplified N-body 
integrations which model the effects of the disc through 
the addition of te rms to cause migration and eccentricit y 
damping (see eg. JSnellgr ove. Papaloizou fc Nelson! i200lf) : 
iLee fc Pealel <2002ft : I Adams fc Lauehlinl <200^ rThe com- 
parison enables us to calibrate these terms and then ap- 
ply the much faster N-body calculations to consider a wider 
range of planet masses and disc surface densities for a longer 
evolution time than is possible for the hydrodynamic simu- 
lations. 

We summarize and discuss our results in the context of po- 
tentially observable commensurabilities in Secti on 6. 

In add ition, we generalize the analytic model of lPapaloizovJ 
(2003) to incorporate migration and circularization effects 
for both planets and also to consider general first order com- 
mensurabilities in an Appendix. Thus this model can be ap- 
plied to systems such as those considered here for which the 
planets have comparable mass and disc interactions may be 
important for either of them. 



2 NUMERICAL SIMULATIONS OF 

MIGRATING PLANETS IN RESONANCE 

We have performed simulations of two interacting planets 
together with an accretion disc with which they also inter- 
act. The simulations perfo rmed here are of the same general 
type a s that performed bv lSnellgrove. Papaloizou fc Nelson! 
( 2001) of the resonant coupling in the GJ876 system in- 
duced by orbital migration caused by interaction with the 
dis c. For details of the numerical scheme and code adopted 
see lNelson et alj ll2000t) . For the simulations performed here 
we adopt n r — 384, and n v — 512 equally spaced grid points 
in the radial and azimuthal directions respectively. 

We use a system of units in which the unit of mass is 
the central mass M» , the unit of distance is the initial semi- 
major axis of the inner planet, T2, and the unit of time is 
27r(GM*/r|) _1,/2 , this being the orbital period on the ini- 
tial orbit of the inner planet. When, as adopted below, r2 
corresponds to 1AU this unit of time is one year. In di- 
mensionless units the inner boundary of the computational 
domain was at r = r m i„ = 0.33, and the outer boundary at 
r — r ma x = 4. In all simulations, the disc model is locally 
isothermal with aspect ratio H/r = 0.05 and the kinematic 
viscosity is set to zero. 



2.1 Initial configuration and computational set up 

Our initial set up shown in Figure ^includes the central star 
with a mass M, and two orbiting planets with masses mi 
and m,2 respectively. The two planets are embedded in the 
disc which is a source of planet orbit migration. They are ini- 
tialized on circular orbits with the central mass taken to have 
a fixed value of 1 solar mass. The gravitational potential was 
softened with softening parameter b = 0.8H. This results in 
the formation of an equilibrium atmosphere around the em- 
bedded planet which then does not accrete. This softening 
also allows an adequate repres entation of type I migration in 
two dimensional discs (see eg. lNelson fc Papaloizoul ([2004 ) ) . 
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Figure 1. The initial configuration: two planets with masses mi 
and ni2, respectively, in circular orbit around a central star with 
mass M* at distances 7*1 and T2 are embedded in a gaseous disc. 



The disc in which planets are initially embedded has the 
initial surface density, E(r), profile specified by 

0.1E (15(r - r min )/r min + 1) if r > r mi „ and 
E(r) = ^ Eo if r > 8r min /5 and 
E (4.5r mi 



,/r) 1 - 5 if r > 4.5r„ 



where r m i„ is the inner edge of the computational domain. 
The planets are located in the flat part of this distribution. 
We use four different values for the maximum value of the 
surface density, E , namely 2 x 10 3 (5.2AU/r 2 ) 2 kg/m 2 , the 
standard value attributed to the minimum mass solar nebula 
at 5.2AU which we denote Ei, then E . 5 = 0.5Ei, E 2 = 2Ei 
and finally E4 = 4Ei. The radial boundaries were taken to 
be open. 



3 TWO DIMENSIONAL HYDRODYNAMIC 
SIMULATIONS: A SURVEY 

We have performed simulations of two interacting planets 
embedded in a disc with which they also interact. This in- 
teraction leads to spiral wave excitation, energy and angular 
momentum exchange between the planetary orbits and the 
gaseous disc which in turn results in orbital migration and 
eccentricity damping. We have carried out simulations for a 
variety of planet masses, initial orbital separations and ini- 
tial surface density scaling in order to explore the possible 
outcomes. Because these simulations are very computation- 
ally demanding, the extent of the survey is limited. Nonethe- 
less some characteristic features emerge. We describe some 
of these results below. 



3.1 One pair of planets in three different 
resonances 

In order to investigate how the orbital evolution depends 
on the surface density of the disc in which the planets 
are embedded we have considered two planets with masses 
mi = 4M® and 7712 = 1M®. We assume that previous evolu- 
tion brought the two planets to a configuration with circular 
orbits with n/ra = 1.2. This separation is slightly smaller 
than that required for a strict 4:3 commensurability. The 



faster migration expected for the outer planet (see equation 
will make the distance between the planets smaller and 
as a consequence of such convergent migration it is possible 
that they ma y become locked in a first order mean motion 
resonance feg. lGoldreicbl(ll965l^ such as 5:4, 6:5, 7:6, 8:7, 
(p + 1) : p, ... and subsequently migrate together maintain- 
ing the commensurability for a considerable period of time 
or perhaps indefinitely, depending on how close the system 
is to a stochastic regime. 



3.2 Attainment of commensurability 

iPapaloizoul <|2003T) has given a simple approximate analytic 
solution for two migrating planets locked in a 2:1 commen- 
surability. Only the outer planet was presumed to interact 
with the disc. As a result of the resonant interaction, the 
eccentricities of both planets grew with time until the ef- 
fects of circularization due to disc tides balanced this. For 
the situation studied here, both planets are embedded in the 
disc and may have significant interaction with it. Further- 
more, higher p commensurabilities than 2:1 occur for planet 
masses in the Earth mass range. Thus in the appendix we 
generalize the solution so that it applies to general first or- 
der commensurabilities and allows for disc tides to act on 
both planets. In particular we show that when the eccen- 
tricities stop growing, provided they are not too large, they 
must satisfy 



e 2 ^ e 2 , 7712711 ai 

td tc2 77ll7l2Ct2 



ei 



t<:2 



f 



migl 



f 

tmig2 J 3 



(5) 



where / = m20i/((p + l)(w.2ai + miai))- 

Here the semi-major axes and eccentricities of the two 
planets (i — 1,2) are at and e;. The migration rates (as- 
sumed directed inwards) and circularization times induced 
by the disc tides are t m i gi = \rii/hi\ — |2aj/(3di)| = 2r r /3, 
and td = \ei/ii\ respectively. The mean motions are n;. 

In addition near a p + 1 : p resonance, the resonant an- 
gles <j> = (p + l)Ai — p\ 2 - H7i, V = (p + l)^i - P^2 - ro 2 
and w\ — zu2 librate about equilibrium values which could 
be near to or n mod 27r. Dissipative effects may be re- 
sponsible for some shift. Here the mean longitudes and the 
longitudes of pericentre of the two planets (i — 1,2) are 
Ai and zui respectively. For giant planets like GJ876 these 
all librate about values close to zero. However, for lower 
mass planets the libratio n of mi — zu2 in particular may 
be ab o ut % mod 2% (eg. ISnellgrove. Papaloizou fc Nelson! 
j200ll) : lLee fc Peald i2002ft t. All the simulations carried out 
here indicate libration about a value closer to 7r than 0. 
Which resonance is established depends on the rate of rel- 
ative migration. Roughly speaking one expects that locking 
occurs only if the relative migration is slow enough that the 
time to migrate through the resonance is longer that a char- 
acteristic libration period. The resonances become stronger 
as p increases, so that increasing the relative migration rate 
tends to cause the attainment of higher p commensurabil- 
ities. However, if the p becomes high enough (p ~ 8 for 
the planet masses considered) the planets enter a chaotic 
region (see eauations l3l4H ) of phase space making commen- 
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5000 io+ 1.5x10* 2x10* 

time (years) 

Figure 2. The evolution of the ratio of semi-major axes for the 
two planets with masses mi = 4Af^ and m,2 = lMg. Starting 
from the lower curve and going upwards the curves correspond to 
the initial surface density scalings So = £4, So = S2, So = Si, 
and So = So. 5 respectively. 



surabilities unstable in the long term (although a system 
may remain in the vicinity of one for a considerable time) 
and introducing sensitivity of detailed outcomes to initial 
conditions. In such cases a scattering may ultimately occur 
that interchanges the positions of the planets. Thus very 
high surface densities with very rapid migration rates do 
not favour attainment of very stable commensurabilities. 
In order to study these effects, we have considered the pair of 
planets evolving in discs with Eo = So. 5, Eo = Ei, £0 = E2 
and Eo = £4. The evolution was followed until a resonance 
was established. The results are summarized in Figure [5] 
where the evolution of the ratio of the semi-major axes, 
which starts from the value 1.2 for all four cases, is shown. 
The fastest migration (steepest slope) corresponds to the 
case where the two planets with mi = 4M^ and 7712 = M9 
are embedded in a disc with Eo = E4 and the slowest to 
a case with a disc with Eo = E0.5. The planets in the disc 
with E = E4 become trapped in 8:7 resonance. If the disc 
surface density is two times smaller then a 7:6 resonance is 
attained, if it is four or eight times smaller then the reso- 
nance attained is 5:4. These results are fully consistent with 
the idea that higher p resonances are associated with faster 
relative migration rates. 

However, this situation corresponds to the evolution dur- 
ing first 18000 years and it is not obvious that the planets 
will remain in the same resonances in their subsequent evo- 
lution. This is because in addition to possibly being in a 
chaotic regime, the ratio of migration time to local orbital 
period varies with disc parameters and accordingly disc lo- 
cation. Consequently, the relative migration might become 
relatively faster as the evolution proceeds resulting in a shift 
to a higher p resonance. We comment that a slowing of rela- 
tive migration, as long as it does not lead to a reversal, will 
not result in a transition to a lower p resonance. Therefore 



changes in effective local migration rates will tend to lead 
to higher p commensurabilities and ultimately chaos. 
The whole evolution for the fastest migration rate is illus- 
trated in Figure El We show there the evolution of the indi- 
vidual semi-major axes, eccentricities, angle between apsidal 
lines and one of the resonant angles. A surface density con- 
tour plot of the disc near the end of the simulation is given 
in Figure |21 together with a comparison between the initial 
and final surface density profiles. 

The inner planet did not migrate significantly until 8:7 res- 
onance was attained at about 6000 years. Subsequently the 
two planets migrated inwards together maintaining the com- 
mensurability. The system appears to be close to the chaotic 
regime as estimated by use of equations j^l-Q). During 18000 
years of evolution the outer planet changed its location from 
a dimensionless radius of 1.2 to 0.9. The eccentricity of the 
outer planet increased slightly and that of the inner planet 
substantially, reaching at around 10000 years the balanced 
value of 0.04 and at later times oscillating around this value. 
At the end of simulation shown here the ratio of eccentrici- 
ties, ei/e2, is equal to 0.125. The local peaks in the values of 
inner planet eccentricity occurring at 2000, 4000 and around 
5400 years correspond to the planet passing through the 5:4, 
6:5 and 7:6 resonances respectively. After about 8000 years 
the angle between the apsidal lines is oscillating around 212° 
and the resonant angle <f> around 264° . The amplitude of the 
oscillations for both angles is decreasing with time. Similar 
behaviour can be seen in the evolution of these planets when 
embedded in a disc with lower surface density (see Figure 

n. 

When Eo = E2, leading to a slower rate of migration, as 
before, the inner planet did not migrate significantly until a 
commensurability was achieved and maintained. A 7:6 reso- 
nance was reached after 10000 years and subsequently main- 
tained. The passage through the 5:4 and 6:5 commensura- 
bilities are clearly marked by a local increase of inner planet 
eccentricity and also in behaviour of the angle between apsi- 
dal lines. In fact after the passage through 5:4 resonance the 
eccentricities do not completely decay potentially making 
the outcome of subsequent resonanc e passages probabilistic 
llKarv. Lissauer fe Greenzweidll99&t) . In this case the plan- 
ets relative migration is slower and the planets stay longer 
in the vicinity of these resonances. At the end of the evolu- 
tion the eccentricities ratio is ei/e2 =0.14. The eccentricity 
of the inner planet is slightly higher than it was in the case 
with Eo = E4. However, the eccentricities continue to grow 
and the equilibrium is not quite established at the end of 
the simulation. The amplitude of oscillations in both the 
angle between the apsidal lines and the resonant angle <j>, is 
smaller than in the previous case. For the two simulations 
that started with even lower disc surface densities such that 
Eo = Ei and Eo = E0.5, the evolution is shown in Figures 
IBTdI 

When Eo — Ei, the evolution of two planets leads to lock- 
ing in 5:4 resonance. This happens at around 8000 years. 
Before that the inner planet migrated very slowly. After be- 
coming trapped in the commensurability the inner planet 
semi-major axis changes in a characteristic oscillatory man- 
ner. These oscillations are clearly visible in all other related 
plots (eg. Figure |SJ . The eccentricity ratio at the end of 
the simulations is ei/e2 = 0.18. The eccentricities are still 
growing. Similar behaviour is seen in angle between the ap- 
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Figure 3. The evolution of the semi-major axes, eccentricities, angle between apsidal lines and the resonant angle for the two planets 
with masses, mi = 4M(j and rri2 = Mj migrating towards a central star, embedded in a disc with initial surface density scaling 
So = £4, obtained by hydrodynamical simulations (four upper panels). The mean surface density profile of the disc near the end of 
the simulation (solid line) and the initial surface density profile (dashed line) and a surface density contour plot are given in two lower 
panels. 
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sidal lines, which starts to librate at time 8000 years around 
180° and at time 18000 years oscillates around 212°. The 
resonant angle changes have very large amplitude and seem 
to grow with time, which might indicate that the planets 
will not remain in this resonance. In this context note that 
equation Q indicates that t m i gi ni increases inwards for a 
uniform surface density. This should favour resonance lock- 
ing at smaller radii. 

Also the planets in a disc with Eo = £0.5 en d in a 5:4 com- 
mensurability. It occurs after approximately 14000 years. 
Because of slow migration rates in low surface density discs, 
the subsequent evolution could not be continued for long 
enough to be able to make a definitive statement about the 
final outcome. However, the oscillations in the resonance an- 
gles are much smaller and more closely centred around 7r 
than is the case when Eo = Si. This could indicate greater 
stability of the 5:4 resonance in this case. 



3.3 Comparison with the analytic model 

It is of interest to compare the eccentricities obtained above 
when the planets are trapped in a commensurability with 
what is expected when resonant effects and disc tides are in 
balance. In the simple analytic model given in the appendix, 
when the eccentricities are small these satisfy equation Jjyl. 
In principle orbital circularization acting on both planets 
should be included. However, for the simulations presented 
here, e-z » ei and even though the circularization time may 
be smaller for the possibly more massive outer planet, the 
effect can be neglected in comparison to that associated with 
the inner planet. Accordingly we set ei = 0. Then equation 
gives the simple relation 



e 2 tm-xn\a\ 
t C 2 \mm2a2 



trriigl tmig2 



(6) 



where we recall that / = m2ai/((p + l)(m2ai + 777,102)). 
We further simplify matters by assuming that p is large. 
Then we may set a\ = 02, ni = 77,2 and neglect / on the 
left hand side of ©. In the same spirit, we replace it by 
m 2/((p + l)(m 2 + rni)) on the right hand side. Then we 
obtain 



1 



mi 



tmigl tmig2 I 3(p + l)(m2 + mi) ' 



(7) 



Interestingly equation JJ) shows that it is the rate of con- 
vergent migration of the two planets that determines the 
eccentricities which are small when this is small in agree- 
ment with our results. Further, other things being equal, 
the eccentricities decrease with increasing p. We also note 
that use of equations Q and gives t c = (T r /W c )(H/r) 2 . 
Thus we obtain after recalling that in general the e folding 
rate for mean motion is a factor of 1.5 greater than that for 
radius, or t mig — 2r r /3 



2 

e 2 = 



m 1 

7712 



-V /2 -i 

a.2 



mi 



0.578(p + l)(m 2 +mi) 



(H/r) 2 .(8) 



We apply this to the case illustrated in Figure [3] for the two 
planets with masses, mi = 4M ffi and 7712 = M ffi embedded in 
a disc with initial surface density scaling Eo = E4 and obtain 
e2 = 0.037 in reasonable agreement with the simulations. 
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Figure 7. The semi- major axis ratio when mi = 4Af®, 7712 = 
lMgj with Eq = So. 5 - upper curve and So = S4 - lower curve. 



3.4 Dependence on the initial separation of the 
planets 

Here we investigate the dependence of attained commen- 
surabilities on the initial radial separation of the planets. 
To do this we have performed simulations with planets of 
the same mass as above but starting with orbital separa- 
tions a little larger than that required for strict 3:2 and 
4:3 commensurabilities. In the former case the two planets 
were initiated with n = 1.32 and r 2 = 1.00. Two initial 
surface density scalings were considered, namely Eo = E0.5 
and Eo = E4. The evolution of the ratio of the semi-major 
axis ratios for both cases is plotted in Figure|7| For Eo = E4 
the end state has the planets in a 9:8 resonance. The evo- 
lution for Eo = E0.5 is shown for comparison. The result 
in the Eo = E4 case is a clear indication that the system 
gets into the chaotic regime as a value of p of the attained 
resonance is increased when the initial semi-major axis ra- 
tio was larger. Stability would have implied that the same 
resonance should have been attained. Unfortunately the mi- 
gration rate in the Eo = E0.5 case was too slow for us to be 
able to attain a commensurability with the available com- 
putational resources and confirm or otherwise the stability 
of the 4:3 commensurability. 

When the ratio of the initial semi-major axes is changed 
from 1.2 to 1.32 the two planets embedded in a disc with 
Eo = E4 become trapped in a different resonance (compare 
Figure |3| and Figure [HJ. This is again consistent with the 
presence of stochastic behaviour as a different migration his- 
tory leads to different end states and it suggests that having 
crossed different resonances before attaining a current one 
may influence the final outcome. 

A similar experiment has been performed but now placing 
the planets close to 4:3 resonance. In this case we have cho- 
sen discs with lower surface density, such that Eo = E0.5 
and So = Ei. In the former case the system enters into 4:3 
resonance while in the latter the system passes through the 
4:3 resonance and becomes trapped in the 6:5 resonance. 
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Figure 9. The semi-major axis ratio for mi = 4Mg and 7712 = 
lMj with surface density scaling parameter So = So. 5 - upper 
curve and So = Si - lower curve. 



Figure 12. The semi-major axis ratio when mi = 7712 = 4Mj) 
with So = Si (uppermost curve) and with So = S4 (two lower 
curves) 



When Eo = E0.5 the planets enter 4:3 commensurability at 
a time around 6000 years. The eccentricities of both planets 
grow. The resonant angle oscillates around 260° with large 
amplitude, higher than it was in the case when the initial 
planet separation was smaller, namely 1.2 (Figure 0. 
For the disc with the surface density Eo = Ei we can follow 
in detail the passage of the planets through the 4:3 resonance 
and the temporary trapping in the 5:4 resonance. In Figure 
[TP at time 4000 years there is a characteristic rapid increase 
in the eccentricities of both planets followed by a slower de- 
crease. The angle between the apsidal lines also shows the 
expected behaviour during resonance crossing. This evolu- 
tion should be compared with that shown in Figure^] for the 
same pair of planets embedded in the same disc but starting 
with a smaller initial planet separation. The effect of a prior 
resonance passage before approaching and becoming rela- 
tively stably trapped in the 5:4 resonance is not present in 
this case. So it is likely that the mutual interaction of planets 
during resonance crossing influences their subsequent evolu- 
tion and whether higher p resonances are attained or not and 
that these planets are either close to or in the chaotic regime. 
This behaviour may be related to the fact that eccentricities 
excited during resonance passages do not completely decay 
between them. 



3.5 Slower relative migration - planets with equal 
masses 

As described above, the resonant interaction must balance 
the tendency towards relative migration of the two plan- 
ets. This is smaller when the planets have the same mass. 
Thus lower p commensurabilities and less tendency to be 
driven into a chaotic state might be expected in this case 
when compared to the situation where the inner planet has 
significantly smaller mass. In order to investigate this, two 
planets of mass 4Me were initiated close to 3:2 resonance 



with ai = 1.32 and 02 = 1 in a disc with Eo = Ei. The 
evolution of the semi-major axis ratio for the planet orbits 
is shown in Figure 1121 We also performed simulations for 
the same masses starting with ai = 1.23 and ai = 1.00 in a 
disc with Eo = E4 and ai = 1.2 and 02 = 1 in a disc with 
Eo = E4. The planets become trapped in the nearest avail- 
able resonance which is a good indication of stability. The 
evolution of the equal mass pair of planets is shown in Fig- 
ures 1131 Tol The planets attaining 3:2 resonance after around 
4000 years show the following behaviour. Their eccentrici- 
ties increased until at 8000 years e\ — 0.013 and e2 = 0.007. 
They then start to oscillate. At around time 13000 the inner 
planet eccentricity dropped to zero for short period of time. 
The resonant angle oscillates around 200°. 
The planets placed close to 4:3 resonance become locked 
in this resonance at around 2000 years. The eccentrici- 
ties of both planets change in a similar way. They do not 
grow much, towards the end of the simulation they oscil- 
late around 0.008. The angle between apsidal lines and the 
resonant angle oscillate around 180° and 230°, respectively. 
The last experiment of this series was as follows. The two 
planets are initially located just interior to 4:3 resonance. 
After 12000 years they arrive through convergent migration 
at the 5:4 resonance and become trapped in it. The eccen- 
tricities of the two planets evolve together in an oscillatory 
way. At time around 16000 years they are very close to zero. 
Particular features at that time are seen in the evolution of 
the angle between apsidal lines as well as the resonant angle 
which are not properly denned for zero eccentricity. 

We also performed calculations with two planets of 30 
M® starting with a\ = 1.3 and 0,2 = 1. The planets were 
embedded in a disc with Eo = Ei. They attained a stable 
3:2 resonance. The behaviour of the ratio of semi-major axes 
close to the resonance is shown in Figure [Tol and other as- 
pects of the evolution are shown in Figure tTTI The excursions 
around resonance are lager and noisier in this case. Note too 
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Figure 11. As for Figure IS1 but when the planets are initiated with a±/a2 = 1.23. Note that the system fails to become trapped in the 
4:3 and 5:4 commensurabilities but nonetheless shows associated strong resonant interaction for long periods of time around t = 4000 
and t = 14000 years respectively. 
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Figure 13. The evolution of semi-major axes, eccentricities, angle between apsidal lines and resonant angle for two planets with equal 
mass, mi = m,2 = 4Mj, migrating towards a central star and embedded in a disc with So = Si. (four upper panels). The mean surface 
density profile of the disc near the end of the simulations (solid line), together with the initial surface density profile (dashed line) and 
a surface density contour plot near the end of the simulation are given in the lowest left and right panels respectively. 
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Figure 17. The evolution of semi-major axes, eccentricities, angle between apsidal lines and resonant angle for two planets with 
masses, mi = 30M^ and m,2 = 30M^ migrating towards a central star embedded in a disc with So = Si (four upper panels). The 
mean surface density profile of the disc near the end of the simulation (solid line) together with the initial surface density profile (dashed 
line) and a surface density contour plot near the end of the simulation are given in lowest left and right panels respectively. 
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4.1 Comparison between hydrodynamic 
simulations and N-body calculations 




^ I i i i i I i i i i I i i i i I i i i i 

5000 10* 1.5x10* 2x10* 

time (years) 

Figure 16. The ratio of semi-major axes for two equal mass 
planets of 30M^, and disc surface density scaling Eq = Sj. 



that in this case the disc planet interaction becomes non 
linear with significant perturbation to the underlying disc 
surface density. There is a tendency towards gap formation 
and the excavation of disc material into two mounds at the 
gap edges. The inner mound occupies only a small radial re- 
gion and so may not be well represented with the available 
numerical resolution in this case. 



4 N-BODY INVESTIGATIONS - A SURVEY 

In order to study migrating planets with a wider range 
of masses and disc surface densities, we have performed a 
resonance survey using a n N-body code. The approach is 
the sa me as that used by [Sji^llgrove^F^ap aloizou fc Nelson! 
l|200lf) and iNelson fc Papaloizoul ll2 02h and it has also 
been u sed bv lLee fc Peald j2002fr and lKlev. Peitz fc Brvdenl 
( 2004). The reader is referred to those papers for details. The 
procedure is to model the effects of a disc by incorporating 
orbital migration and eccentricity damping as described by 
equations 10- through the addition of appropriate accel- 
eration terms to the equations of motion. We summarize the 
parameters and outcomes of some of the N body simulations 
in Table □ 

We begin by discussing a calibration procedure which en- 
ables a matching of the N body simulations with the hydro- 
dynamic simulations. The advantage of the N body simula- 
tions is that the time evolution of a system can be followed 
for a much longer period of time. After discussing a few 
examples, we move on to discuss the results, shown in Ta- 
ble Q on the possible commensurabilities to be expected for 
low mass planets for particular choices of their masses, disc 
parameters and migration and eccentricity damping rates. 



The hydrodynamical calculations discussed in the previous 
Section allow us to follow planets migrations for almost 
2xl0 4 years. As a result of this we have been able to simu- 
late planets becoming trapped in resonances. The behaviour 
and stability of the resonance trapping varies with the planet 
masses and the surface density of the disc in which they are 
embedded. The outcome of these simulations could be well 
matched to those of N-body integrations where we incorpo- 
rate the simple prescriptions for the migration and eccen- 
tricity damping given through equations 11121 . 
Using these expressions in the N-body code we have ex- 
tended the hydrodynamic calculations for a longer period of 
time and studied the long term stability of the resonances 
we have found in that context. As a first step we have ad- 
justed the numerical coefficients in equations 11121 in such a 
way that the hydrodynamic and N-body approach give the 
same qualitative evolution. 

As an example in Figure ITH1 we show the results of the com- 
parison for the case of two planets with masses 1 Mq and 
4 respectively embedded in a disc with initial surface 
density scaling parameter So — £4. The numerical coeffi- 
cients adopted were W m = 0.3647 and W c = 0.225. These 
results show good agreement with the hydrodynamic results 
(see Figure |3J and display the same 8:7 commensurability. 
The fitted coefficients were also reasonably close to those ex- 
pected from the analytic disc planet interaction theory, con- 
firming both the migration and eccentricity damping times. 
Though the latter were typically ~ 40% longer than ex- 
pected from the analytic theory. The evolution was followed 
using the N-body approach for an additional 1.2 x 10 6 years. 
The two planets remain in the 8:7 resonance during this time 
and there is no indication of any significant changes in the 
monitored quantities for the last 9 x 10 5 years. The long 
term evolution of this system is shown in Figure 1191 It is 
interesting to note that the equilibrium value of eccentricity 
of the inner planet is around = 0.03, in accordance with 
the value predicted by the simple analytic model discussed 
in Section 2 and derived in the Appendix. The angle be- 
tween apsidal lines is close to 180° the other resonant angles 
occupy a wide band between 80° and 287°. A similar pro- 
cedure has been applied to other cases presented in Section 
3. The values of the numerical coefficients (W m ,W c ) used 
in equations (11121 obtained by comparing and matching the 
planet evolution calculated by hydrodynamic and N-body 
codes are given in Tabled 

In the case of two planets with masses, mi = 4Mj and 
m.2 = lMj migrating towards a central star embedded in 
the disc with Eo = £2 the situation is less stable. The 7:6 
resonance found in the hydrodynamic simulations is main- 
tained for around 2 x 10 s years, later there is a shift into a 
10:9 commensurability. Finally at around 6 x 10 years scat- 
tering occurred. Also the 5:4 commensurability in the lower 
surface density case (Eo = Ei) was maintained only until 
1.2 x 10 s years, then there is a transition to 6:5 which lasted 
till 2.1 x 10 5 years. Next, the planets evolved till 6 x 10 5 years 
locked in 7:6 resonance and after that continued in 8:7. That 
was the status of the evolution at the end of our calculations 
at 1.3 x 10 6 years. The same pair of planets embedded in 
the disc with even lower surface density (Eo = E0.5) mi- 
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Table 1. This table lists some of N body calculations of interacting planets with imposed migration and 
eccentricity damping rates, assumed to result from interaction with a disc, that we performed, together 
with their outcomes. During these calculations the migration has been followed for a period of time in 
which the inner planet changes its semi-major axis from «2 = 1 to (12 = 0.385 in dimcnsionlcss units. 
The outer planet is initially placed in circular orbit at r = 1.4 in dimensionless units. Thus the situation 
corresponds with appropriate scaling to the inner planet being initiated at 5.2AU and then migrating to 
2AU. Column 1 gives the planet mass ratio, columns 2 and 3 give the outer and inner planet mass in 
Earth masses respectively. Columns 4-7 indicate the outcome of the interaction as either attainment of a 
commensurability or a scattering for those cases where the migration was too rapid for that to occur. Where 
there is a scattering, the radius where it occurred is indicated in dimensionless units in brackets. An asterisk 
indicates that the resonance showed signs of instability that would result in the transition to a higher p 
resonance on further inward migration which could not be followed. 
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Table 2. The fitting parameters obtained by comparing hydrodynamic and N-body calculations 
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Figure 18. The evolution of semi-major axes, eccentricities, angle between apsidal lines and resonant angle for two planets with 
masses, mi = 4Mgj and mi = IMq migrating towards a central star embedded in the disc with S = S4 obtained by N-body simulation 
(dotted lines in two upper panels and dark lines in two lower panels) superimposed onto the hydrodynamic simulation shown already 
in Figure 131 



grated together in 5:4 resonance till 1.15 x 10 6 years, then 
they moved into a 6:5 commensurability. These outcomes are 
suggestive o f the long term instability of the high p commen- 
surabilities (iKarv. Lissauer fc Greenzweidll993tl . 

We have noticed a strong sensitivity of outcomes to 
the fitting parameters (especially in the numerical coeffi- 
cient Wm- This is symptomatic of chaotic motion. A very 
small change in the values of the numerical coefficients given 
in Table |21 can produced large variations in outcomes. The 
value W c is difficult to fit for pairs of two equal mass planets 
because of smaller and less regular changes in eccentricities 
obtained in the hydrodynamic simulations of such systems 
(see for example Figure Hoi . 



4.2 N-body resonance survey 

We applied the much faster N-body calculations to extend 
our hydrodynamic simulations to a wider range of planet 
masses and disc surface densities. The results are sum- 
marized in Table We allowed the numerical coefficients 



Wm and W c in equations IUI2H . when used in the N body 
evolution, to be free parameters that we could choose to 
provide the best match between the hydrodynamic and N 
body methods. In t his way we could test the applicabil- 
ity of the result s of j Tanaka. Takeuchi fc Ward! i2002t) and 
iTanaka fc Ward! feOuj) to this problem. 

When we chose W m = 0.3704 in equation Q and W c = 
0.289 in equation for the case when the planets have 
equal masses, (m-i/mi) = 1, their initial separation in cir- 
cular orbits in dimensionless units is such that rxjri = 1.4 
and they are embedded in a disc which extends from r = 0.3 
to r — 1.4, than the most likely outcome of their orbital 
evolution caused by the disc-planet interaction, in agree- 
ment with the hydrodynamic simulations, is attainment of 
a 3:2 commensurability. This is the case for a wide range of 
disc surface density scalings. For all values of Eo considered 
here, this commensurability was attained. A transition from 
the 3:2 commensurability to the 4:3 commensurability for a 
given mass ratio (bigger than 1) is expected to occur when 
either the surface density scaling and/or the mass ratio is 
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Figure 19. The evolution of the semi-major axes, eccentricities, angle between apsidal lines and resonant angle for two planets with 
masses, m\ = 4Mg) and ni2 = lMg migrating towards a central star embedded in the disc with £ = £4 obtained with N-body 
simulations. 



large enough. The value of p for the resonance which might 
be established accordingly increases. 

When mi/7712 > 6 scattering is likely to occur rather than 
stable resonance trapping, especially for higher surface den- 
sities. This is clearly indicated in Table Q 



5 CONCLUSIONS 

We have investigated migration induced resonant capture in 
a system of two planets for the most part in the mass range 
1 — 4A% embedded in a gaseous disc. We also considered a 
system with equal masses of 30M®. As the disc was laminar, 
apart from in the 30M^ case, the planets underwent type I 
migration which occurs through the linear response of the 
disc to perturbation by the planets. We considered disc sur- 
face densities in the range 0.5 — 4 times that expected for 
the minimum mass solar nebula at 5.2 AU. Thus migration 
rates should thus be typical of those expected for standard 
protoplanetary discs. 

Our conclusions are based on the results of two surveys in 
which we studied the evolution of a pair of planets with a 



range of planetary masses (1 - 30M®), disc surface densities 
and initial planet separations. The smaller survey used hy- 
drodynamic simulations and because of computational re- 
quirements was limited to run lengths of ~ 2 x 10 4 inner 
planet orbits. Similarly initial conditions which started with 
the planets too far apart were precluded. A more extensive 
survey used a simplified N-body approach where the influ- 
ence of the disc was modelled through the addition of exter- 
nal torques and eccentricity damping terms (see equations 
and ©). 

The numerical calculations are supplemented with an ana- 
lytic model describing two planets migrating in resonance 
with arbitrary eccentricities. We found a very good agree- 
ment between numerical and analytical results for an estab- 
lished commensurability. 

Both our hydrodynamic and N-body approaches give clear 
evidence of an interesting relation between the mass ratio in 
a planetary pair and the type of the resonance which is ex- 
pected to be established as a consequence of tidally induced 
orbital migration. When both planets have similar masses, 
provided the disc surface density at their locations is very 
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similar, the rate of convergent or relative migration is small. 
This favours attainment of low p commensurabilities. An 
example of this occurred in the equal 4M® case we consid- 
ered. There, the planets became locked in the nearest first 
order commensurability lying between them. Thus if they 
were assembled with 1.55 > ri/ra > 1.31 the migration led 
to trapping in 3:2 resonance. 

Here, we comment that the very well known example of 
such a system, namely the two largest mass planets orbit- 
ing a millisecond radio pulsar PSR B1257+12, are close to 
a 3:2 commensurability. The most recent determination of 
the masses of these plane ts gives them as 4.3 ± 0.2M m and 
3.9±0.2M© respectively, jKonacki fc Wolszczanll2003ft . The 
mass ratio is thus very close to unity and according to our 
findings if migration in a standard quiescent protoplanetary 
gaseous disc is responsible for the evolution then a possi- 
ble outcome could be attainment of a 3:2 commensurabil- 
ity. The orbital p eriods for both planets are 66.5419 days 
and 98.2114 days jKonacki fc Wolszczanll2003l) which places 
them near 3:2 resonance. Such a resonance could be formed 
and maintained through the processes discussed in this pa- 
per, with the planets later moving slightly out of resonance. 
More specific detailed modeling will be presented elsewhere. 
For more disparate mass ratios, our simulations indicate the 
attainment of first order commensurabilities with higher p 
(for example 4:3, 5:4, 7:6 and 8:7). The dynamics develops 
a stochastic character, and these may be maintained for a 
different periods of time before becoming unstable. Then of- 
ten a passage from one resonance to another takes place. 
The short duration of the hydrodynamic simulations pre- 
cluded extensive exploration of stochastic behaviour or the 
long term stability of resonances. However, we could investi- 
gate these issues employing the simplified N-body approach. 
In this context we comment on a result we have obtained 
for an 8:7 commensurability, being the stable end state of 
the evolution of two planets in circular orbit with the initial 
separation ratio ri/r 2 = 1.2, and mass ratio m\/m 2 = 4 
and the highest disc surface density scaling considered here 
with So = £4. The stable configuration was maintained 
until 1.4 x 10 6 orbits at which point our calculations were 
stopped. Of course this result should be considered in the 
context of extreme sensitivity to initial conditions and other 
parameters. Nonetheless, the final value of the inner planet 
eccentricity is in good agreement with the prediction of our 
analytic model for established commensurabilities. 
Thus while stochastic behaviour tends to disrupt higher p 
commensurabilities, the calculations presented here do not 
enable us to rule them out entirely. 

Simple N-body calculations allowed us to extend our inves- 
tigations not only to consider longer term behaviour but 
also a wider range of initial conditions. These calculations 
strengthen the conclusions derived from the hydrodynamic 
simulations. 

Finally, we comment further on the onset of stochastic be- 
haviour which we encountered both in our hydrodynamic 
simulations and N-body calculations. The reason for this 
phenomenon can be related to resonance overlap. Accord- 
ing to our estimates we might expect isolated resonances in 
which a system of planets in the mass range considered can 
be locked and migrate together to be such that the integer 
p in the expression p + 1 : p for a first order commensura- 
bility is < 8. For larger p, as we can see from Table a 



system is likely to undergo a scattering and exchange of or- 
bits of the two planets. However, stochastic behaviour may 
extend to values as small as 4. Stable commensurabilities 
with sufficiently small p correspond to physically possible 
planetary configurations and they should be the centre of 
further investigation. 

Among all planned future space missions COROT will be 
the first with the capability to observe planets in the mass 
range relevant to our investigation. Studies of commensu- 
rabilities have already been applied successfully to analyse 
the motion of the pulsar planets, and they proved to be a 
powerful tool in that context. They have the potential to be 
similarly useful in our search for Earth-like planets in other 
systems. Detection of such resonances can also yield useful 
information about orbital migration as a process operating 
during planet formation. 
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APPENDIX 

A simple model 

The equations of motion for a system consisting of 2 planets 
and primary star, modeled as point masses moving under 
their mutual gravitational attraction, are conveniently ex- 
pressed inl^amiltonian form using Jacobi coordinates (see 
ISinclairl jl975l) ). Then the radius vector r 2 , of the inner 
planet of reduced mass m 2 is measured from the primary 
star of mass M 4 and that of the outer planet, n, of reduced 
mass mi is referred to the centre of mass of the primary star 
and inner planet. The Hamiltonian can be written correct 
to second order in the planetary masses as 



H 



1, ,. ,2 , ,. ,2x GM.imi GM, 2 m 2 

-(mi n + m 2 r 2 J : : : 1 — 

2 |ri| |r 2 

Gm\m 2 Gm\m%r\ ■ r 2 



(9) 

|ri 2 | |n |" 

Here M*i = M* + mi,M» 2 = Af* + m 2 and n 2 = r 2 - n,. 
The Hamiltonian can be expressed in terms of the osculat- 
ing semi-major axes, eccentricities and longitudes of perias- 
tron a,i,ei,zui,i = 1,2 respectively as well as the longitudes 
Ai, and the time t. We recall that Xi — m(t — ten) + Wi, 
with rti = yj GM,i /a? being the mean motion and tot giv- 
ing the time of periastron passage. The energy is given 
by Ei = —GmiM*i/(2ai), and the angular momentum 
hi — rriiyj GM„iai(l — ef) which may be used to describe 
the motion instead of o» and e;. 

The Hamiltonian may quite generally be expanded in a 
Fourier series involving linear combinations of the three 
angular differences X j — zu i , i = 1,2 and vji — vj 2 (eg. 
iBrouwer fc Clemencd (119611) ). 
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Near a first order p + 1 : p resonance, we expect that 
both <p — (p+ l)Ai — p\2 — tui, and ip = (p + l)Ai — p\2 — ra?2, 
will be slowly varying. When resonances are non overlap- 
ping so that the p + 1 : p resonance may be isolated, in its 
neighborhood, terms in the Fourier expansion involving only 
linear combinations of <j> and ip as argument are expected 
to produce the largest perturbations, while the remainder, 
which will have rapidly varying argument, may be averaged 
out. 

The resulting Hamiltonian (oc mimz) may then be written 
in the general form 



Hi 



a i 



C k ,i(ai/a 2 , ei, e 2 ) cos(fc0 + lip), (10) 



where in the above and similar summations below, the sum 
ranges over all positive and negative integers (fc, I) and the 
dimensionless coefficients Ct,i depend on ei, e 2 and the ratio 
0,1/0,2 only. Here we shall not need to specify these further. 
We also make the inconsequential simplification of replacing 
M»i by M m . 



Basic Equations 

We take the Hamiltonian system derived from (1101 and mod- 
ify it by incorporating torques and rates of change of energy 
that act on each of the protoplanets and which are presumed 
to be derived from interaction with the disc. 

The equa t ions o f motion are very similar to those given 
in lPapaloizoul J2003) . They differ only in that the additional 
forces derived from disc interaction are allowed to act on 
both protoplanets rather than only the outer one in this 
case. 

dEi/dt = -mdHi 2 /d\i - (mTi/y/1 - ef + A), 
dhi/dt = -dH 12 /d\i - dHii/dwi - Ti, 
dXi/dt = rn + mdHvi/dEi + dHw/dhi, 
d-uji/dt = 8Hi2/dhi, 

The external disc torque acting on the planet mi is — Ti. 
Associa ted wi th this, we remove orbital energy at a rate 
ruTi I yl — ef from irn which would correspond to the action 
of a disc density perturbat ion rotating with a pattern speed 
nil \f \ — e? (see also ISnellgrove. Papaloizou fc Nelson! 
l|200l|h iNelson fc Papaloizoul (120021) '). We include an addi- 
tional energy loss rate Di for rrii which leads to an orbital 
circularization time for mi given by 



tci — 



mmiejy/ GM»o» 
A(l-e?) ■ 



(11) 



This circularization time is such that if the the other planet 
was absent such that rrn was affected only by the central 
star and disc, et would decay according to 



dd 
~dt 



(12) 



We presume that Ti, and Di are given as functions of the 
orbital parameters aj, ej associated with m 4 . In addition we 
suppose that they lead to orbital evolution on a timescale 
sufficiently long compared to other effects that they may 
always be considered to act as small perturbations to the 
equations of motion. However, we do not assume an expan- 
sion in terms of small eccentricities. 



dm 
~dT 



We thus obtain to lowest order in the perturbing masses. 
3(P+lWm 3 ^ Cfc<i(jfe + 0sinW + iV>) 



Af» 

3niai 
GM*mi 



(13) 



dn2 
~~dt 



3pn 2 mia2 
M7ai 

3n 2 (2 2 

GM,m 2 



^C fc ,i(fc + I) sin(k</> + lip) 



ri2T 2 



+ D 2 



(14) 



dei 
~dt 



d 7712711 

T7i _ 



eiM„ 



C Ki sin(fc0 + li>) (k-{p+ l)(fc + 0(1 - Vl-e?)) (15) 



dt2 _ e 2 7nia 2 n 2 -\/r 



dt t C 2 aie 2 M* 

^2c kil Bin(k<p + lip) (l+p(k + l)(l- v/l - e|)) (16) 

^ = (p + l)m - pn 2 - ^2(D ktl + E k ,i) cos(fc0 + lip). (17) 

^ = (p+ l)m -pn 2 -^(D fc ,i + F M ) ms(kcp + lip). (18) 
Here 

2(p + l)nia?m 2 9 
T>kA = ~ — — (yh,i/ai) 



M* dai 



2pn 2 a 2 mi d 
Ml da~2~ 



(C h ,,/ai) ; 



(19) 



E k 1 



and 



Fk.i 



nirri2((p + 1)(1 - ef) - Py /l - ef ) dC fc „ 



eiM* 



9ei 



pn 2 a 2 mi(-y/l — — 1 + e 2 ) 9C fc! ; 



aie 2 M* 



<9e 2 



(20) 



(p + l)mm 2 (l -ej-y/l- ej) dC k ,i 
eiM, Sei 
: 772a 2 mi((p +1) ^/l - e 2 - p (l - ej)) gg^j 



aie 2 M* 



<9e 2 



We further note that <p— ip — Z02 — roi is the angle between 
the two apsidal lines of the two planetary orbits. 



Time Independent Solutions 

When no external disc torques or dissipation act (T; = Di — 
0) time independent or stationary solutions may occur for 
which ip, (p and each of ni,n 2 ,ei,e 2 are constant. In prin- 
ciple any stationary values of ip = ipo and <p = <po might be 
possible. Obvious possibilities, as also indicated in our nu- 
merical simulations, are values of <po and ipo that are either 
0, or 7r. In that case it follows directly from equations (1131161 
with T — Di = that there are stationary solutions with 
77i,n2, ei, e 2 being constant. In other cases, consideration of 
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equations 1131 - 1181 indicates that we should regard 0o and 
ipo as functions of ei,e2 and the ratio 02/01. 
A relation between ei,e2 and the ratio 02/01 for steady so- 
lutions follows by subtracting equations 1171 and (11811 in the 
form 

cos(fc0 o + li>o) = ^ cos(fc0 o + JV>o) = S p .(22) 

This condition in fact matches the precession rates of the 
orbits of the two planets such that they maintain a fixed 
orientation. It then additionally follows from equations l|17[l 
and itTHl that 



(p + 1)77.1 - pn 2 = D k,i cos(fc0 o + lipo) + S p . 



(23) 



This gives a further relation between ei,e2 and the ratio 
aa/oi. We could for example specify ei and then both ei and 
0,2/0,1 would be specified but in any case the quantity ((p + 
l) n i)/(P n 2) — 1, supposing mi and 777,2 to be comparable, is 
at least of order of smallness of the mass ratio mi/M, <C 1. 
Thus in a perturbation quantity that is already first order in 
the perturbing masses we may set 02/01 = ((p + l)/p) _2//3 
which gives the condition for a p + 1 : p commensurability. 



Time Dependent Solutions with Migration and 
Circularization 

We now generalize the above solution to incorporate the ef- 
fects of small non zero disc torques Ti and dissipation Di. 
In this case the two planets migrate inwards locked in res- 
onance with n\/n2 maintained nearly equal to p/(p + 1). 
In the absence of any tidal effects which could act to cir- 
cularize the orbits the eccentricities are found to increase 
monotonically with time. We note that equations 1131141 
define inward migration timescales or e folding rates for m 
as t m igi = GM,m,i^Jl — / (STiOiUi) . We suppose that in 
order to accommodate small values of Ti and Di , the angles 
tpi = tp — tpo and <f>i — — 0o, which measure the departures 
from the values appropriate to the steady state solution, 
take on non zero values of small magnitude which can also 
change very slowly with time. These magnitudes are pre- 
sumed sufficiently small that we may employ a first order 
Taylor expansion so that the perturbation to sin(fc0 + lip) is 
(fc0i + Itpi) cos(fc0o + lipo) and similarly the perturbation to 
cos(fc0 + lip) is — (fc0i + Itpi) sin(fc0o + ■ 
When the migration time is very long compared to any other 
time scale in the problem, equations 1171 and 1181 are of 
the same form as in the time independent case, apart from 
small corrections proportional to the migration rate which 
we neglect. These then lead to the same conclusions as in the 
time independent case. Accordingly we conclude that there 
are relationships between ei , e2 and a 2 /a\ which specify the 
latter two once e\ is specified but that always 02 /a\ satisfies 
the condition for (p + 1) : p resonance with a correction of 
order of the perturbing masses. 

Equations Q13I16JI then become 



dni 3(p + l)nlm2 

~dT = Ml ^ 

3niai niTi 
GM^rni ^/T~e 



+ 



C M (fc + Z)(fc0i +lipi) 



+ Dl 



(24) 



dn 2 
~~dT 



?>P"n\m\a2 
M t ai 



^C w (Hi)Wi+if) 



3n20 2 
GM* 7712 



W2T2 



0" 



+ D 2 



(25) 



dei 
~dt 



ei 



m 2 n ly Jl — e\ 



eiM„ 



X ^C M (fe0i + lipi) (k-(p+ l)(k + 0(1 - Vi-e?)) ( 26 ) 



de 2 _ e 2 niia 2 n 21 /l - e 
~dt ~ 



t,::2 



aie 2 M, 



xJ^CkAWi + lif>i) (l+p(k + l)(l- V 1 - 6 !)) • (27) 

Here Ck,i = Ck,i cos(fc0o + Zt/>o). 
Taking the ratio of 1241 and 1251 gives 

dni 
dn 2 

( P + l)nfm a aiX)«fc.j(* + + (-7&T + D >) 

rA (28) 

pn\mi02 E n k ,(k + 1)- + D2 j 

with TLk.i = Ck,i{k<f>i + lipi). 

The ratio of 1261 and 1271 gives another relation of the form 



m2e 2 aini 

c? de 2 

E «*,«(' + V{k + 0(1 - y/1 ~ el)) + 



e 2 M * a 1 



-e^mia 2 n 2 t c 2 



E«k,j(* - (p + i)(fc + 0(1 - V 13 ^?)) + /tA ; 

Using pri2 = (p + 1)tt-i, equation 1281 gives 

</>i 53 fcc fc ,i(fc + + Vi^ ;Cfc ' i ( fc + ^ = 

-(p + l)ma?m 2 (77=? + °i) + pmm2aia 2 (-^=f + ^2 

Gmiffi2 ((p + l) 2 n\m,2ai + p 2 n\m\a2) 

Given that we may regard e2 and 02/01 as functions of ei, 
we may also regard equations l|28|l and l|29^ as determining 
and tp as functions of ai, and ei. Using this information, 
equations 1241 and 1261 may be used to determine n\ and 
ei as functions of time so completing the solution. 
One readily finds that these satisfy 

1 dni 



(29) 



(30) 



m dt 



miGt2 



(it^T + + aim2 + (^=fc 



[m 2 ai + mio 2 ] 



(31) 



dei 



ei 



((37^7 + 



(l-ef)tel 



eie2(A — 1) dei 



At c2 de 2 
) ~ (377^ + 



+ 



T 



A [msoi + 771102] 



(32) 



where 
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ei de 2 
ei dei 



(33) 



p + 1 



1 - (p+ 1)(1 - A /l-e?)-p A / 1 - e 2+P) (34) 



A = 1 + 

and 
T 



Equation I3H implies that ni always increases with time 
corresponding to inward migration. Notably this equation 
does not depend on the order of the resonance. 
Equation 132H governs the evolution of the eccentricities. It 
indicates that when e\ = e-z = 0, given that the derivative 
de 2 /dei remains finite, the eccentricities will increase with 
time provided that l/t m i g \ > l/tmig2- This is just the con- 
dition for the migration of the uncoupled two planet system 
to be convergent, in the sense that the relative separation 
should decrease. When the circularization times are finite, 
there is the possibility of steady state eccentricities found 
by equating the right hand side of (1321 to zero. In the limit 
that they are small, the equilibrium eccentricities must thus 
satisfy: 



e| e 2 m2ni ai 

t c i t c 2 77li77,2a2 



2 2 

r-r) / = 

Icl Ec2 



1 



migl 



mig2 



f 



(35) 



where / = m,2a 1 /((p + l)(m2tli + mia 2 )). Note that when 
the migration and circularization times t m i g and t c scale in 
the same way as the mean motions at an orbital location, 
the equilibrium eccentricities are independent of orbital lo- 
cation. 
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